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Master curves for the stress tensor invariants in stationary states of 
static granular beds. Implications for the thermodynamic phase space 

Luis A. Pugnaloni,^ Jose Damas,^^ Iker Zuriguel,-* Diego MazaP 

We prepare static granular beds under gravity in different stationary states by tapping the 
system with pulsed excitations of controlled amplitude and duration. The macroscopic 
state — defined by the ensemble of static configurations explored by the system tap after 
tap — for a given tap intensity and duration is studied in terms of volume, V, and force 
moment tensor, E. In a previous paper [Pugnaloni et al., Rhys. Rev. E 82, 050301(R) 
(2010)], we reported evidence supporting that such macroscopic states cannot be fully 
described by using only V or E, apart from the number of particles A'^. In this work, 
we present an analysis of the fluctuations of these variables that indicates that V and 
E may be sufficient to define the macroscopic states. Moreover, we show that only one 
of the invariants of E is necessary, since each component of E falls onto a master curve 
when plotted as a function of Tr(E). This implies that these granular assemblies have a 
common shape for the stress tensor, even though it does not correspond to the hydrostatic 
type. Although most results are obtained by molecular dynamics simulations, we present 
supporting experimental results. 



I. Introduction 

The study of static granular systems is of funda- 
mental importance to the industry to improve the 
storage of such materials in bulk, as well as to op- 
timize the packaging design. However, far from 
yielding such benefits of practical interest, physi- 
cists have found a fascinating challenge on their 
way that has stuck them, to a large extent, in the 
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area of static granular matter. Such challenge is 
finding an appropriate description of the simplest 
state in which matter can be found: equilibrium 
[l]. At equilibrium, a sample will explore differ- 
ent microscopic configurations over time, in such a 
way that macroscopic averages over large periods 
will be well defined, with no aging. Moreover, if 
not only equilibrium but ergodicity is present, av- 
erages over a large number of replicas at a given 
time should give equivalent results to time aver- 
ages [2]. Finally, one expects that all macroscopic 
properties of such equilibrium states can be put in 
terms of a few independent macroscopic variables. 
Then, a thermodynamic description and, hopefully, 
a statistical mechanics approach can be attempted. 
Theoretical formalisms based on these assumptions 
have been used to analyze data from experiments 
and numerical models. However, some of the foun- 
dations are still supported by little evidence. 
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The use of careful protocols to make a granular 
sample explore microscopic configurations within a 
seemingly equilibrium macroscopic state has given 
us the first standpoint [3||^. In such protocols, the 
sample is subjected to external excitations of the 
form of pulses. Well defined, reproducible time av- 
erages are found after a transient if the same pulse 
shape and pulse intensity are applied. However, for 
low intensity pulses, a previous annealing might be 
in order since equilibrium is hard to reach — as in 
supercooled liquids. We will use the expressions 
equilibrium state, steady state or simply state to 
refer to the collection of all configurations gener- 
ated by using a given external excitation after any 
transient has disappeared. 



More than twenty years ago, Edwards and Oak- 
shott [rj put forward the idea that the number of 
grains N and the volume V are the basic state vari- 
ables that suffice to characterize a static sample of 
hard grains in equilibrium. The NV granular en- 
semble was then introduced as a collection of mi- 
crostates, where the sample is in mechanical equi- 
librium, compatible with N and V. However, newer 



II. Simulation 



theoretical works [8 -[13 suggest that the force mo- 
ment tensor, E, (S = Va, where a is the stress ten- 
sor) must be added to the set of extensive macro- 
scopic variables (i.e., an NVT, ensemble) to ade- 
quately describe a packing of real grains. 



In the rest of this paper, we show experimen- 
tal and simulation evidence that the equilibrium 
states of static granular packings cannot be only 
described by V (or equivalently the packing frac- 
tion, (j), defined as the fraction of the space covered 
by the grains) nor by S. We do this by generat- 
ing states of equal V but different E and states of 
equal E but different (j). We also show that states 
of equal V may present different volume fluctua- 
tions. Moreover, we show that states of equal V 
and E display the same fluctuations of these vari- 
ables, suggesting that no other extensive parameter 
might be required to characterize the state (apart 
from TV). Finally, but of major signiflcance, we 
show that the shape of the force moment tensor 
is universal, in the sense that different states that 
present the same trace of the tensor actually have 
the same value in all the components of E. 



We use soft-particle 2D molecular dynamics 14|15 



Particle-particle interactions are controlled by the 
particle-particle overlap ^ = d — \rij\ and the ve- 



locities r,; 



and ojj. 



Here, represents the 



center-to-center vector between particles i and j, d 
is the particle diameter and a; is the particle an- 
gular velocity. These forces are introduced in the 
Newton's translational and rotational equations of 
motion and then numerically integrated by a ve- 
locity Verlet algorithm 16 . The interaction of the 



particles with the fiat surfaces of the container is 
calculated as the interaction with a disk of infinite 
radius. 

The contact interactions involve a normal force 
F„ and a tangential force F^. 



Fn = fcnC - InV. 



Ft = -min(/i|F„|,|F,|)-sign(C) 



where 



(t') dt' 



2 J 



(1) 



(2) 



(3) 



(4) 



(5) 



The first term in Eq. ([T]) corresponds to a restor- 
ing force proportional to the superposition ^ of the 
interacting disks and the stiffness constant The 
second term accounts for the dissipation of energy 
during the contact and is proportional to the nor- 
mal component w"^- of the relative velocity r^j of 
the disks. 

Equation ([2| provides the magnitude of the force 
in the tangential direction. It implements the 
Coulomb's criterion with an effective friction fol- 
lowing a rule that selects between static or dynamic 
friction. Notice that Eq. ^ implies that the max- 
imum static friction force |-Fs| used corresponds to 
//|f„|, which effectively sets ^idynamic = ^static = A^- 
The static friction force Fs [see Eq. (|3|] has an 
elastic term proportional to the relative shear dis- 
placement C and a dissipative term proportional to 
the tangential component vl j of the relative veloc- 
ity. In Eq. (15]) , s is a unit vector normal to . The 
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elastic and dissipative contributions are character- 
ized by kg and 7s, respectively. The shear displace- 
ment C is calculated through Eq. Q by integrating 
vjj from the beginning of the contact (i.e., t = to). 
The tangential interaction behaves like a damped 
spring which is formed whenever two grains come 
into contact and is removed when the contact fin- 
ishes |17j . 

The particular set of parameters used for the 
simulation is (unless otherwise stated): fi — 0.5, 
kn = 10^{mg/d), 7„ — 300{m^/gjd), kg = and 
7s — 200(m^/gjd). In some cases, we have varied 
fi and 7„ in order to control the friction and resti- 
tution coefficient. The integration time step is set 
to 6 ~ 10^'^ y^d/g. The confining box (13.39(i-wide 
and infinitely high) contains N — 512 monosized 
disks. Units are reduced with the diameter of the 
disks, d, the disk mass, to, and the acceleration of 
gravity, g. 

Tapping is simulated by moving the confining 
box in the vertical direction following a half sine 
wave trajectory [Asin(27rj/t)(l — 0(27rzyi— tt))]. The 
excitation can be controlled through the amplitude, 
A, and the frequency, v, of the sinusoidal trajectory. 
We implement a robust criterion based on the sta- 
bility of particle contacts to decide when the sys- 
tem has reached mechanical equilibrium 1 14 before 
a new tap is applied to the sample. Averages were 
taken over 100 taps in the steady state and over 20 
independent simulations for each value of A and v. 

The volume, V, of the system after each tap 
can be obtained from the packing fraction, 0, as 
V = NTr{d/2)'^ /(f). We measure in a rectangu- 
lar window centered in the center of mass of the 
packing. The measuring region covers 90% of the 
height of the granular bed (which is of about 40c?) 
and avoids the area close to the walls by 1.5d. We 
have observed that 4> is sensitive to the chosen win- 
dow. However, none of the conclusions drawn in 
this paper are affected by this choice. 

The stress tensor, a, is calculated from the 
particle-particle contact forces as 

where the sum runs over all contacts. 

The force moment tensor, S, is defined as S = 
Va. During the course of a tap S is non-symmetric, 
however, once mechanical equilibrium is reached in 



accordance with our criterion, S becomes symmet- 
ric within a very small error compared with the 
fiuctuations of S. Although E may depend on 
depth, we have measured the force moment tensor 
by simply summing over particle-particle contacts 
in the entire system. 

The fluctuations of cj) (Acj)) and S (AE) are cal- 
culated as the standard deviation in the 100 taps 
obtained in each steady state. We average cj), E, 
and their fluctuations over 20 independent runs for 
each steady state and estimate error bars as the 
standard deviation over these 20 runs. 

III. Experimental method 




Figure 1: (a) Schematic diagram of the experimen- 
tal set-up. A: accelerometer, S: shaker, C: camera, 
O: oscilloscope, FG: function generator, PC: com- 
puter, (b) Example of an image of the packing. 
The region of measurement is indicated by a red 
square. 

The experimental set-up is sketched in Fig. [ija). 
A quasi 2D Plexiglass cell (28 mm wide and 150 mm 
high) is filled with 900 alumina oxide beads of di- 
ameter d = 1 ± 0.005 mm. The separation between 
the front and rear plates was made 15% larger than 
the bead diameter. The cell is tapped by an elec- 
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tromagnetic shaker (Tiravib 52100) with a trend 
of half sine wave pulses separated five seconds be- 
tween them. The tapping amplitude was controlled 
adjusting the intensity, A, and frequency, i/, of the 
pulse and was measured by an accelerometer at- 
tached to the base of the cell. Averages are taken 
over 500 taps after equilibration. 

High resolution images (more than 10 MPix) are 
taken after each tap. To calculate the packing frac- 
tion, we only consider a rectangular zone at the 
center of the packing whose limits are at 4 mm 
from the borders [see Fig. [l|^b)]. We determine the 
centroid of each particle by means of a numerical 
algorithm with subpixel resolution. Then, we cal- 
culate the packing fraction by assuming that the 
2D projection of each bead corresponds to a disk 
of diameter d = 1 mm. We estimate the packing 
fraction with a resolution of ±0.001. As the sep- 
aration between plates is larger than the particle 
diameter, small overlaps between the spheres are 
present in the 2D projections of the images. There- 
fore, the calculated packing fraction in some dense 
configurations might result slightly higher than the 
hexagonal disk packing limit (7r-\/3/6). 



actual pulse used to drive the system is merely a 
control parameter but not a macroscopic variable 
that describes the state. Therefore, the external 
pulse does not need to be described with a simpli- 
fied quantity. The complete functional form of the 
pulse can be given instead. In our case, we use a 
sine pulse and both, the pulse amplitude and fre- 
quency, are needed to fully describe the excitation. 
We will employ the usual nondimensional peak ac- 
celeration F = Upeak/g — A{2iTvY / g (whcrc g is 
the acceleration of gravity) and the frequency u to 
precisely define the external excitation. A detailed 
study of the dynamical response of a granular bed 
to a pulse of controlled intensity and duration can 
be found in Damas et al. 



22 



One major issue in studying equilibrium states 
is the evidence that one can have, indicating if the 
system is actually at equilibrium [5, . Since the defi- 
nition of equilibrium is circular |23| , we can simply 
do our best to check if different properties of the 
system have well defined means (as well as higher 
order moments of the distributions) which should 
not depend on the history of the processes applied 
to the sample. 



Tapping characterization and 
asymptotic equilibrium states 



IV. 



It is often debated [18j[19| what is the appropri- 
ate parameter to characterize the external excita- 
tion used to drive a granular sample. Dijksman et 
al. 18 proposed a parameter related to the lift- 
ofT velocity of the granular bed. Ludewing et al. 



19 presented an energy based parameter. Pug- 



naloni et al. 15 20 suggested that the factor of 



expansion induced on the sample would be a suit- 
able measure [21]. The usual perspective to define 
a pulse parameter, e, is to achieve a collapse of the 
<t)-e curves as different details of the pulse shape 
are changed (such as amplitude and duration) . Pa- 
rameters defined in all these previous works fail to 
make such curves collapse for the data presented 
in this paper. The main reason for this is that our 
4>-e curves are non-monotonic (see for example Fig. 
|6]) presenting a minimum whose depth depends on 
the details of the pulse shape. Therefore, a sim- 
ple rescale of the horizontal axis does not suffice to 
collapse the curves. 

Since we are interested in macroscopic states, the 
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Figure 2: Evolution of towards the steady state 
starting from a disordered configuration initially 
obtained by strong taps (experimental results). Re- 
sults for two tapping intensities are shown: F — 4.8 
(blue), and F = 17 (green). In both experiments, 
the frequency of the pulse is v — 2>Q Hz. 

We have generated configurations corresponding 
to a particular pulse (of a given shape and inten- 
sity) by repeatedly applying such pulse to the sys- 
tem and allowing enough time for any transient to 
fade. To prove that our samples are at equilib- 
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rium, we approach a particular pulse through differ- 
ent paths — and starting from different initial con- 
ditions (ordered and disordered configurations) — 
and confirm that the steady states obtained present 
equivalent mean values and second moments of the 
distributions of the variables of interest. 
In Ref. 



0.92 
0.9 
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we showed that the steady states 
corresponding to excitations of high intensity are 
reached in a few taps, even if the initial configu- 
ration corresponds to a very ordered structure. In 
Fig. [2] we consider the equilibration process from a 
highly disordered initial configuration. For low tap- 
ping intensities (blue line in Fig. [2]), the packing 
fraction evolves to the steady state in two stages. 
Just after switching the tapping amplitude to the 
reference value, the system rapidly evolves to values 
of 4> close to the final steady state. Beyond this ini- 
tial convergence, a slower compaction phase takes 
the system to the final steady state. For high tap- 
ping intensities, the evolution to the steady state is 
very rapid. The steady state is reached after about 
a hundred taps [green line in Fig. [2])]. Therefore, 
we apply a sequence of at least 1000 taps in all our 
experiments before taking averages to warrant that 
the steady state has been reached. In our simula- 
tions, 400 taps of equilibration were enough. 

An interesting example of equilibration is pre- 
sented in Fig. |3] In this case, we present the values 
of (j) and of Tr(I]) during a very special sequence 
of pulses obtained in our simulations. The system 
is initially deposited from a dilute random config- 
uration with all particles in the air. First, we tap 
the system 200 times at F = 4.9, then 800 times 
at F = 61.5 and finally, 200 times at F = 4.9. In 
the whole run, we keep v = 0.5y/ g/d. These two 
values of the pulse intensity have been chosen be- 
cause they are known to produce packings with the 
• in the steady state 



same mean (p m tue steady state 1^. However, the 
mean E is clearly different. Notice, however, that 
the same values of 0, E and its fluctuations, A(f) 
and AE, are observed for the steady state of low F 
obtained before and after the 800 pulses of high F. 

Unless otherwise stated, all the results we present 
in what follows correspond to steady states. We 
have tested this by obtaining the same states 
through two different preparation protocols consist- 
ing of: (i) application of a large number of identi- 
cal pulses starting from a disordered configuration, 
and (ii) application of a reduced number of iden- 
tical pulses after annealing the system from much 
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Figure 3: Evolution of cj) (a) and Tr(E) (b) as the 
pulse intensity is suddenly changed between two 
values that produce steady states with the same 
(j) but different E in our simulations. The initial 
200 taps correspond to F = 4.9, the middle 800 
pulses to F = 61.5 and the final 200 pulses, again, 
to F = 4.9. In all cases, v — Q.b^Jg/d. The mean 
values and standard deviations in each section are 
indicated by arrows. 



higher pulse intensities. In a few cases, the results 
(mean values and/or fluctuations) from both pro- 
tocols did not match. This was an indication that 
a steady state, if existed, was not reached by one 
of the protocols (or by both protocols). This hap- 
pens especially for low intensity taps, which require 
longer equilibration times. Such cases have been 
removed from the analysis. 

V. Volume and volume fluctuations 

In Fig.jljja), we plot in the steady state as a func- 
tion of F from our experiments for different pulse 
frequencies [24] . As we can see, there exists a min- 
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imum (j) at relatively high T. A similar experiment 
on a three-dimensional cell also yielded analogous 
results 26 . This behavior has also been reported 

An explanation for 



15 20 25 



for various models 

this, based on the formation of arches 
given in [15 



v(Hz) 



(a) 



has been 

The position of the minimum in (f) 
shifts to larger F if the frequency of the pulse is 
increased (i.e., if the pulse duration is reduced). 
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Figure 4: (a) Experimental results for the steady 
state packing fraction, 0, as a function of the re- 
duced peak acceleration, F, for different frequen- 
cies, f, of the tap pulse, (b) Histogram of the con- 
figurations visited by the system for = 30 Hz: 
F = 15 (green) and F = 28 (red). 

Although the increment in the packing fraction 
beyond the minimum is rather small, it is impor- 
tant to remark that this difference is not an artifact 
introduced by our experimental resolution. In Fig. 
|4|^b), we show the histogram for the sequence of 
packings obtained for F ~ 15 (the minimum pack- 
ing fraction for ly — 30 Hz) and for F = 28 (the 
largest excitation explored for = 30 Hz). Both 
steady states are statistically comparable; however, 
it is possible to distinguish different mean values. 

Since steady states of equal cj) obtained at both 




Figure 5: The steady state volume fraction fluctua- 
tions, A0, as a function of F from our experiments 
with v — 30 Hz. The solid line is only a guide to 
the eyes. 



sides of the </> minimum are generated via very dif- 
ferent tap intensities, it is worth assessing if such 
states are, in fact, equivalent. This can be done by 
comparing the volume fluctuations of such states. 
A similar analysis was done in Ref. 



27 for states 



generated with diff'erent pulse amplitude and dura- 
tion in liquid fluidized beds. 

The fluctuations of (j) in the steady state, as mea- 
sured by the standard deviation A</), are presented 
in Fig. [5] As we can see, a minimum in the fluctu- 
ations is just apparent. The position of such min- 
imum in the fluctuations coincides with the mini- 
mum in (p. Unfortunately, the resolution of (j) in our 
experiments is around the size of the fluctuations. 
However, the results from the simulations are not 
limited in this respect and we address to those for 
a more reliable assessment of the fluctuations. 

In Fig. ^a) , we plot in the steady state as a 
function of F for our simulations. Although the 
fluctuations of (j) are large [see Fig. |3], its mean 
value is well defined with a small confidence inter- 
val (see error bars). For low excitations, 4> decreases 
as F is increased. However, beyond a certain value 
Fmin, the packing fraction grows. The same trend 
is observed if the tap frequency v is changed. How- 
ever, the minimum is deeper for lower v and its 
position F,„in shifts to larger values of F as is 
increased in coincidence with our experiments (see 
Fig. |4| . Due to the change in the depth of the cj) 
minimum in the (jj-T curves, a simple rescaling of F 
is unable to collapse the data for different frequen- 



cies. However, rescaling (j) and F with 



and 
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^min, respectively, yields very good collapse, even 
between simulated and experimental data 



26 



In Fig. |6f^b), the volume fraction fluctuations, 
A(/), are plotted as a function of F. As we can see, 
these fluctuations are non-monotonic, as suggested 
by our experiments (Fig. [s]). Non-monotonic vol- 
ume fluctuations have also been reported in Ref. [6] . 
For A(j), we obtain a minimum and a maximum. We 
have also observed a maximum in for values of 
r below the ones reported here (see for example 
Refs. f26l and [25]). However, we do not report 
such low values of F in this work and focus on tap- 
ping intensities that warrant the steady state with 
a modest number of pulses. 
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Figure 6: (a) Simulation results for the volume frac- 
tion, ^, as a function of the reduced peak accelera- 
tion, F, for different frequencies, of the tap pulse, 
(b) The corresponding volume fraction fluctuation, 
A(/), as a function of F. 

The value of F, at which the fluctuations dis- 
play a minimum, coincides with Vmim the value at 
which the minimum packing fraction, (^mim is ob- 
tained. The maximum coincides with the inflection 
point in the 0-F curve at higher F. Since one ex- 
pects to find few mechanically stable configurations 



compatible with a large volume (low (/)), it seems 
reasonable that fluctuations reach a minimum if </) 
does so. Similarly, there should be few low-volume, 
mechanical stable configurations which implies that 
at high fluctuations should diminish. Hence, a 
maximum in should be present at intermedi- 
ate packing fractions. This is more clearly seen in 
Fig. [T] where we plot the fluctuations as a function 
of the average value of 0. 
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Figure 7: Density fluctuations as a function of 
for different frequencies of the tap pulse in the sim- 
ulations. 

Figure [7] presents two distinct branches: the 
lower branch corresponds to F > F„jj„ and the up- 
per branch to F < F.,„i„. For F > Fmi„, the fluc- 
tuations corresponding to different tap durations 
collapse, suggesting that such equal-(/), equal-A0 
states might correspond to unique states (below, we 
will find that this is not the case) . For F < Fmin , we 
obtain states of same as some states in the lower 
branch but presenting larger fluctuations. This is 
clear evidence that the equal-0 states of the upper 
and lower branch are indeed distinct and that other 
macroscopic variables must be used to distinguish 
one from the other. 

We have assessed a number of other structural 
descriptors (coordination number, bond order pa- 
rameter, radial distribution function). In all cases, 
equal-0 states from the upper and lower branch 
of the 0-A()f) curve (Fig. [?]) present similar values 
of the structural descriptors with only subtle dis- 
crepancies. Although this indicates that the states 
are not equivalent, it also suggests that such de- 
scriptors are not good candidates to form a set of 
macroscopic variables, along with V ^ to uniquely 
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identify a given steady state. 

In Ref. [26], we have assessed the force moment 
tensor S as a good candidate to complete V and N 
in describing a stationary state. This has been sug- 
gested by some theoretical speculations |12 . How- 
ever, some authors prefer to directly replace V by 
E. Below, we will show that both, V and E, are 
required for the equilibrium states generated in our 
simulations. Moreover, we will show that only one 
of the invariants of E is necessary (at least in our 
2D systems) and that fluctuations of these variables 
suggest that no other extra macroscopic parameter 
may be required. 

VI. Stress tensor 

Before we focus on the force moment tensor, we 
will consider the stress tensor, cr, in order to under- 
stand the phenomenology of the force distribution 
in our tapped granular beds. We recall that E and 
a are simply related through E = Vcr. However, 
we have to bear in mind that V is not a simple 
constant since the volume of the system depends, 
in a nontrivial way, on the shape and intensity of 
the excitation. 

In Fig. [8j we show the components of u as a func- 
tion of F for different v. As a reference, we show 
results of our simulations for a frictionless system. 
In a frictionless system, the shear vanishes and ayy 
is only determined by the weight of the sample since 
the Janssen's effect is not present. As we can see, 
the frictionless sample presents a constant value of 
ayy for all F. For low F, the frictional samples dis- 
play values of ayy below the frictionless reference. 
This is a consequence of the Janssen's effect since 
part of the weight of the sample is supported by 
the wall friction. Consequently, in this region, axy 
is also positive [see Fig. [8](b)]. However, for each 
V, there is a critical value of F, Tshear=o- Beyond 
it, the sample presents an apparent weight above 
the weight of the packing. In correspondence with 
this, a^y changes sign and becomes negative. This 
indicates that, for F > Tshear=o^ the frictional walls 
are not supporting any weight. Rather, they pre- 
vent the packing from expansion by a downward 
frictional force. As F is increased beyond Tshear=o, 
the packing tends to store most of its stress in the 
horizontal direction (axx) while ayy eventually sat- 
urates. For very intense pulses, the sample expands 



and lifts off significantly during the tap. When the 
bed falls back, it creates a very compressed struc- 
ture with most of the stress transmitted in the lat- 
eral directions and the wall friction sustaining the 
system downwards. It is worth mentioning that 
^shear=o is always higher than Tmin (see Fig. [6|. 
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Figure 8: (a) Diagonal components of the stress 
tensor, a, as a function of F for different frequen- 
cies V of the tap pulse (the upper set of curves cor- 
responds to ayy and the lower set to axx)- (b) Off 
diagonal component of the stress tensor, axy The 
horizontal line corresponds to Oyy [in panel (a)] and 
to axy [in panel (b)] from simulations of a packing 
of frictionless disks. 



VII. The force moment tensor mas- 
ter curve 

Since the stress is not an extensive parameter, the 
force moment tensor is generally used to character- 
ize the macroscopic state 12|13 . Therefore, we will 



use E in the rest of the paper. Let us simply remark 
that since V presents a non-monotonic response to 
F, the curves in Fig. [8] present a somewhat differ- 
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ent shape if E is plotted instead of a. Particularly, 
Hyy does not display a minimum at low F, as the 
one observed for ayy in frictional disks but a mono- 
tonic increase. In Fig. [9] we show the trace of S as 
a function of F. There is a clear monotonic increase 
of Tr(E) as F is increased. Moreover, for a given F, 
if the frequency of the excitation pulse is increased, 
a significant reduction in the force moment tensor 
is observed. 




Figure 9: Trace, Tr(I]), of the force moment tensor 
as a function of F for different frequencies ly of the 
tap pulse. 



In Fig. 10 



we plot the components of S as a func- 
tion of its trace for all the steady states generated. 
We can see that all data for different F and ly col- 
lapse into three master curves. This indicates that 
if two equilibrium states present the same Tr(I]), 
all the components of E are also equal. Here, we 
point out a relevant piece of information that will 
be discussed in the next section. Two states may 
present equal force moment tensor but differ in vol- 
ume. This means that many points collapsing in 
Fig. [T0| correspond to states of different (j>. There- 
fore, at equilibrium, irrespective of the structure of 
the sample, two states with the same trace in E will 
present equal E. 

In a liquid at equilibrium, the stress tensor is 
diagonal and all elements along the diagonal are 
equal. This hydrostatic property allows us to know 
the full stress tensor if we only know the hydro- 
static pressure (i.e., if we only know the trace of 
the tensor). In our granular samples, the force mo- 
ment tensor can also be known if the trace is known. 
However, the shape of the tensor in static packings 



under gravity is defined by the three master curves 
of Fig. [10 
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Figure 10: (a) Diagonal components of the force 
moment tensor, E, as a function of Tr(E) for dif- 
ferent frequencies of the tap pulse (the upper set of 
curves corresponds to Tiyy and the lower to Y^xx)- 
(b) Off diagonal component of the force moment 
tensor, E^y. 

To our knowledge, there is no previous specula- 
tion that this property must hold for static granular 
packings. A more detailed study on the extent of 
this commonality of the shape of the force moment 
tensor will be pursued in a future paper 



28 . How- 



ever, we show some suggestive preliminary results 
below. 

In Fig. [TT] we show the components of E as a 
function of Tr(E) for a range of samples of differ- 
ent materials, for different tapping intensities, and 
for different tapping frequencies. As we can see, 
there is reasonable collapse of the data onto the 
same three master curves shown in Fig. 10 This is 



an indication that these master curves may be uni- 
versal and enclose a rather fundamental underlying 
property (inaccessible to us at this point) of static 
granular beds. 
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Figure 11: (a) Diagonal components of force mo- 
ment tensor, S, as a function of Tr(I]) for different 
frequencies of the tap pulse, different friction co- 
efficient and different restitution (the upper set of 
curves corresponds to and the lower to S^a;)- 
(b) Off diagonal component of the force moment 
tensor, The dashed lines are only a guide to 

the eyes. 



why the force moment fluctuations should present 
a minimum. Provided that the minimum of ATr(S) 
coincides with the minimum of Ac/), it can be spec- 
ulated that a reduced number of geometric con- 
figurations can accommodate a limited number of 
force configurations. We have seen that all individ- 
ual components of E present the same minimum 
in their fluctuations, however, the actual values in 
ATr(E) are dominated by AT^^x which takes values 
five times larger than AS^j,. 
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VIII. Force moment tensor fluctua- 
tions 

Since we have shown in the previous section that 
only the trace of S suffices (along with the master 
curves) to describe the full force moment tensor, 
we will now focus on this invariant and its fluctu- 
ations. In Fig. 12 the fluctuations of Tr(I]) are 



plotted as a function of F. We obtain a single min- 
imum, in contrast with the minimum and maxi- 
mum observed in A^. Interestingly, the states with 
minimum ATr(E) correspond to the states where 
the minimums and At/) are reached for each v. 
However, unlike A(/), the depth of the minimum in 
ATr(E) is fairly independent of v. It is unclear 



Figure 12: (a) Fluctuations of the trace of the force 
moment tensor as a function of F for different fre- 
quencies of the tap pulse, (b) Fluctuations of the 
trace of the force moment tensor as a function of 
Tr(E). 

If we plot ATr(E) in terms of the average value 
Tr(E) [see Fig. [l2](b)], we can see that the curves 
collapse on top of each other over a wide range of 
Tr(E). However, some deviations are apparent at 
very low and very high forces. Although the fair 
collapse of the curves suggests that states of equal- 
E may correspond to the same equilibrium states, 
we will see in the next section that many of these 
equilibrium states are distinguishable through the 
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volume. The inflection observed at high Tr(S) cor- 
responds to the change of regime observed in Fig. 
[Sf^a) where the vertical stress saturates and most of 
the contact forces are directed in the x-direction. 



tions even if they correspond to an equivalent mean 
volume or an equivalent mean force moment tensor 
(see states joined by solid lines in Fig. 14). 



IX. The thermodynamic phase 
space 

As we have suggested, the mean volume of static 
granular samples is not sufficient to describe the 
equilibrium state since states of equal-iji may 
present distinct fluctuations. On the other hand, 
the force moment tensor seems to be able to serve 
as a standalone descriptor since states of equal S 
do generally present the same E fluctuations. How- 
ever, states of equal-S may present different vol- 



umes. In Fig. 13 we plot the loci of the equilibrium 
states generated in our simulations in a hypothet- 
ical (j)-T, thermodynamic phase space. As we can 
see, states of equal V but different S are obtained 
as well as states of equal E and different V. 
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Figure 14: (a) Phase space 0-Tr(E). (b) Fluctu- 
ations of the state variables. Selected neighboring 
states are colored in pairs for comparison. Distant 
states of equal V or Tr(S]) are joined by thick lines. 



Figure 13: Phase space (^-Tr(E). Loci visited in 
the simulations for different frequencies of the tap 
pulse. 

We can ask now if these two state variables suffice 
to fully describe the equilibrium states. A hint that 
this may be the case is given by the fact that states 
generated with different P and v but that corre- 
spond to the same state in the (j)-T, plot display the 
same fluctuations of these variables. In Fig. [Mj we 
have highlighted some pairs of neighboring states. 
We can see that such states do also present similar 
fluctuations [Fig. [T4j[b)]. In contrast, states which 
are distant in the plot present distinct fluctua- 



Let us point out here, that the sole coincidence 
of fluctuations in the macroscopic variables is not 
a rigorous proof that the set of chosen variables 
is a full complete set of thermodynamic param- 
eters. Future explorations of these systems may 
confirm or disprove that N, V and Tr(E) are a 
full set of macroscopic, extensible variables able to 
describe all equilibrium states. Meanwhile, it is 
clear that for moderate tapping intensities, around 
which the minimum in is observed, the approxi- 
mation of a simple NV or N'S ensemble is not war- 
ranted in view of the large discrepancies between 
the curves generated with different pulse frequen- 
cies in Fig. [T3[ 
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X. Concluding remarks 

We have studied steady states of mechanically sta- 
ble granular samples driven by tap like excitations. 
We have varied the external excitation by changing 
both, the pulse amplitude and the pulse duration. 
We have considered the macroscopic extensive vari- 
ables V (volume) and S (force moment tensor) , and 
their fluctuations. From the results, we can draw 
the following conclusions: 

• There seems to be a rather robust set of master 
curves for which implies that the knowl- 
edge of Tr(E) suffices to infer the other com- 
ponents of the force moment tensor. 

• The equilibrium states cannot be only de- 
scribed by y or S, apart from the number of 
particles N. 

• The equilibrium states seem to be well de- 
scribed by the set NVTt{J:). 

There exists a number of points to be considered 
in view of these findings. Here, we mention a few 
that may serve as starting points for future direc- 
tions of research: 

• What is the extent to which the E master 
curves are applicable? Is this dependent on 
the dimensionality of the system, the excita- 
tion procedure, the chosen contact force law, 
etc? 

• What is the dynamics during a single pulse 
that leads to the appearance of the 4> mini- 
mum? Is this minimum present in states gen- 
erated with other types of pulses like fluidiza- 
tion or shear? The ubiquity of this minimum 



in simulation models 15 20 suggests that it 



might be found in numerous conditions. 

• Are the fluctuations shown in Figs. [7] and 
[12] the definitive phenomenological equation of 
states? Other authors have so far found mono- 
tonic density fluctuations 27 or concave up 
density fluctuations [6]. 

• How much of the 0-Tr(I]) plane can be ex- 
plored by changing material properties? 

• Are there other excitation protocols (such as 
shearing) that may give rise to steady states 



that are thermodynamically equivalent to the 
ones obtained by tapping? 

• Is it possible to construct a phenomenological 
entropy function from the equations of states 
(Figs, [t] and [l2| ) by simple integration of a 
Gibbs-Duhem-like equation? Let us bear in 
mind that Fig. [7] is multiply valued. 

It is worth stressing that if two ensembles gen- 
erated by arbitrary excitation protocols — such as 
tapping or shearing — happen to present the same 
mean values (and fluctuation) for all macroscopic 
variables, then such macroscopic states should be 
considered thermodynamically identical. However, 
it may be the case that a given protocol produces 
a narrow range of macroscopic states that can be 
eventually described with a reduced set of macro- 
scopic variables. 
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